In [1]:
#Setup the working environment
setwd("/hpc/group/pbenfeylab/mz187/")
In [2]:
# Load the basic package_Seurat
suppressMessages(library(tidyverse))
suppressMessages(library(Seurat))
suppressMessages(library(cowplot))
suppressMessages(library(ComplexHeatmap))
suppressMessages(library(circlize))
suppressMessages(library(GeneOverlap))
suppressMessages(library(gprofiler2))
suppressMessages(library(ggrepel))
suppressMessages(library(muscat))
suppressMessages(library(limma))
suppressMessages(library(scran))
suppressMessages(library(future))
suppressMessages(library(RColorBrewer))

Warning message:
“package ‘tidyverse’ was built under R version 4.2.3”
Warning message:
“package ‘tibble’ was built under R version 4.2.3”
Warning message:
“package ‘tidyr’ was built under R version 4.2.3”
Warning message:
“package ‘readr’ was built under R version 4.2.3”
Warning message:
“package ‘purrr’ was built under R version 4.2.3”
Warning message:
“package ‘dplyr’ was built under R version 4.2.3”
Warning message:
“package ‘stringr’ was built under R version 4.2.3”
Warning message:
“package ‘forcats’ was built under R version 4.2.3”
Warning message:
“package ‘lubridate’ was built under R version 4.2.3”
Warning message:
“package ‘ComplexHeatmap’ was built under R version 4.2.1”
Warning message:
“package ‘circlize’ was built under R version 4.2.3”
Warning message:
“package ‘GeneOverlap’ was built under R version 4.2.1”
Warning message:
“package ‘gprofiler2’ was built under R version 4.2.3”
Warning message:
“package ‘muscat’ was built under R version 4.2.2”
Warning message:
“package ‘limma’ was built under R version 4.2.3”
Warning message:
“package ‘scran’ was built under R version 4.2.3”
Warning message:
“package ‘SingleCellExperiment’ was built under R version 4.2.2”
Warning message:
“package ‘SummarizedExperiment’ was built under R version 4.2.2”
Warning message:
“package ‘MatrixGenerics’ was built under R version 4.2.1”
Warning message:
“package ‘matrixStats’ was built under R version 4.2.3”
Warning message:
“package ‘GenomicRanges’ was built under R version 4.2.3”
Warning message:
“package ‘BiocGenerics’ was built under R version 4.2.1”
Warning message:
“package ‘S4Vectors’ was built under R version 4.2.3”
Warning message:
“package ‘IRanges’ was built under R version 4.2.3”
Warning message:
“package ‘GenomeInfoDb’ was built under R version 4.2.3”
Warning message:
“package ‘Biobase’ was built under R version 4.2.3”
Warning message:
“package ‘scuttle’ was built under R version 4.2.3”
Warning message:
“package ‘RColorBrewer’ was built under R version 4.2.3”
In [4]:
sessionInfo()

R version 4.2.0 (2022-04-22)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: CentOS Stream 8

Matrix products: default
BLAS/LAPACK: /hpc/group/pbenfeylab/mz187/miniconda3/envs/muscat/lib/libopenblasp-r0.3.21.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] future_1.32.0         ggrepel_0.9.3         gprofiler2_0.2.2     
 [4] GeneOverlap_1.34.0    circlize_0.4.15       ComplexHeatmap_2.14.0
 [7] cowplot_1.1.1         Seurat_3.1.5          lubridate_1.9.2      
[10] forcats_1.0.0         stringr_1.5.0         dplyr_1.1.2          
[13] purrr_1.0.1           readr_2.1.4           tidyr_1.3.0          
[16] tibble_3.2.1          ggplot2_3.4.2         tidyverse_2.0.0      

loaded via a namespace (and not attached):
 [1] Rtsne_0.16          colorspace_2.1-0    rjson_0.2.21       
 [4] ggridges_0.5.4      IRdisplay_1.1       GlobalOptions_0.1.2
 [7] base64enc_0.1-3     clue_0.3-64         leiden_0.4.3       
[10] listenv_0.9.0       fansi_1.0.4         codetools_0.2-19   
[13] splines_4.2.0       doParallel_1.0.17   IRkernel_1.3.2     
[16] jsonlite_1.8.5      ica_1.0-3           cluster_2.1.4      
[19] png_0.1-8           uwot_0.1.14         sctransform_0.3.5  
[22] compiler_4.2.0      httr_1.4.6          Matrix_1.5-4.1     
[25] fastmap_1.1.1       lazyeval_0.2.2      cli_3.6.1          
[28] htmltools_0.5.5     tools_4.2.0         rsvd_1.0.5         
[31] igraph_1.5.0        gtable_0.3.3        glue_1.6.2         
[34] RANN_2.6.1          reshape2_1.4.4      Rcpp_1.0.10        
[37] vctrs_0.6.3         ape_5.7-1           nlme_3.1-157       
[40] iterators_1.0.14    lmtest_0.9-40       globals_0.16.2     
[43] timechange_0.2.0    lifecycle_1.0.3     irlba_2.3.5.1      
[46] gtools_3.9.4        MASS_7.3-60         zoo_1.8-12         
[49] scales_1.2.1        hms_1.1.3           parallel_4.2.0     
[52] RColorBrewer_1.1-3  reticulate_1.25     pbapply_1.7-0      
[55] gridExtra_2.3       stringi_1.7.6       S4Vectors_0.36.0   
[58] foreach_1.5.2       caTools_1.18.2      BiocGenerics_0.44.0
[61] shape_1.4.6         repr_1.1.6          bitops_1.0-7       
[64] rlang_1.1.1         pkgconfig_2.0.3     matrixStats_1.0.0  
[67] evaluate_0.21       lattice_0.21-8      ROCR_1.0-11        
[70] patchwork_1.1.2     htmlwidgets_1.6.2   tidyselect_1.2.0   
[73] parallelly_1.36.0   RcppAnnoy_0.0.20    plyr_1.8.8         
[76] magrittr_2.0.3      R6_2.5.1            gplots_3.1.3       
[79] IRanges_2.32.0      generics_0.1.3      pbdZMQ_0.3-9       
[82] pillar_1.9.0        withr_2.5.0         fitdistrplus_1.1-11
[85] survival_3.5-5      future.apply_1.11.0 tsne_0.1-3.1       
[88] crayon_1.5.2        uuid_1.1-0          KernSmooth_2.23-21 
[91] utf8_1.2.3          plotly_4.10.0       tzdb_0.4.0         
[94] GetoptLong_1.0.5    data.table_1.14.8   digest_0.6.31      
[97] stats4_4.2.0        munsell_0.5.0       viridisLite_0.4.2
In [3]:
# Read rice root atlast Copilot project
rc.integrated <- readRDS("Rice_12S_WT_atlas_20221117_seu3.rds")
In [4]:
table(rc.integrated$orig.ident,rc.integrated$celltype.fine.anno)

        
         Stem Cell Niche Sclerenchyma Root Cap Atrichoblast Trichoblast Cortex
  sc_108             508          194      310          255         266    961
  sc_109             327          176      212          156         785    770
  sc_115             157           92      155          444        1299    373
  sc_116             268          198      240          194        1260   1083
  sc_192             781          765      704          121         123   1791
  sc_193             627          617      585          102         102   1417
  sc_194             404          478      632          120          31   1184
  sc_195             437          509      626          124          44   1209
  sc_196             300          723      245           97         326    933
  sc_7               225          189      134           54          74    911
  tz1               1923         2024      413         2364         485   2626
  tz2               1789         2528      468          644         606   3999
        
         Endodermis Phloem Pericycle Exodermis Protoxylem Metaxylem
  sc_108        558    246       751       381        376       189
  sc_109        690    116       423       435        325        84
  sc_115        136     98       203       243        162        19
  sc_116        210    262       668       274        467        64
  sc_192        250    200       282       656        894        83
  sc_193        152    177       248       496        693        80
  sc_194        144    116       246       537        650        66
  sc_195        143    118       280       536        693        83
  sc_196        123     69       151       976        562        44
  sc_7          175    101       451       184        494        50
  tz1           697    207       671      2782       1625       162
  tz2          1084    745       714      2826        987       364
In [5]:
table(rc.integrated$orig.ident,rc.integrated$timezone.ID.P)

        
         Elongation Maturation1 Maturation2 Meristem
  sc_108        912         465        1339     2279
  sc_109        699        1011         831     1958
  sc_115        374        1279         886      842
  sc_116        879        1523        1043     1743
  sc_192       1805         432        2159     2254
  sc_193       1434         366        1654     1842
  sc_194       1194         298        1635     1481
  sc_195       1224         344        1751     1483
  sc_196       1550         662        1156     1181
  sc_7          713         250         852     1227
  tz1          4531        2296        3729     5423
  tz2          5229        1887        2976     6662
In [6]:
table(rc.integrated$orig.ident)

sc_108 sc_109 sc_115 sc_116 sc_192 sc_193 sc_194 sc_195 sc_196   sc_7    tz1 
  4995   4499   3381   5188   6650   5296   4608   4802   4549   3042  15979 
   tz2 
 16754
In [6]:
# Check the cell type annotation for highly confidently annotated cells
options(repr.plot.width=16, repr.plot.height=16)
DimPlot(rc.integrated, reduction = "umap", group.by = "timezone.ID.P",
        pt.size = 1.0, 
        order = c("Maturation2","Maturation1","Elongation","Meristem"),
        cols = c("#ebf8e3", "#51c8bd", "#009ac8", "#005fa8"),label = F)
In [4]:
zscore <- function(x){(x-mean(x))/sd(x)}
In [10]:
known.good.markers <- read.csv("./1-Root_Atlas/Ethylene_genes.csv")
In [17]:
known.good.markers$Gene.ID
  1. 'LOC-Os07g03250'
  2. 'LOC-Os09g39320'
  3. 'LOC-Os02g32610'
  4. 'LOC-Os04g52140'
  5. 'LOC-Os06g40360'
  6. 'LOC-Os02g10700'
  7. 'LOC-Os07g48630'
  8. 'LOC-Os09g31400'
  9. 'LOC-Os04g38400'
  10. 'LOC-Os07g06130'
  11. 'LOC-Os04g46220'
  12. 'LOC-Os01g58420'
  13. 'LOC-Os08g31580'
  14. 'LOC-Os10g25170'
  15. 'LOC-Os03g08460'
  16. 'LOC-Os03g08470'
  17. 'LOC-Os01g21120'
  18. 'LOC-Os06g09390'
  19. 'LOC-Os09g11460'
  20. 'LOC-Os03g49500'
  21. 'LOC-Os05g06320'
  22. 'LOC-Os04g08740'
  23. 'LOC-Os02g57530'
  24. 'LOC-Os03g12350'
  25. 'LOC-Os07g48430'
  26. 'LOC-Os05g46240'
  27. 'LOC-Os03g58520'
In [7]:
length(known.good.markers$Gene.ID)>1
TRUE
In [8]:
rc.integrated

An object of class Seurat 
96918 features across 79743 samples within 3 assays 
Active assay: integrated (19395 features, 19395 variable features)
 2 other assays present: RNA, SCT
 5 dimensional reductions calculated: pca, umap, umap_50, umap_3D, umap_2D
In [18]:
## Check gene existence in single cell data
match(known.good.markers$Gene.ID,rownames(rc.integrated@assays$SCT@data))
  1. 30364
  2. 29059
  3. 10671
  4. 18000
  5. 22583
  6. 9789
  7. 25422
  8. 28570
  9. 17087
  10. 23539
  11. 17623
  12. 2973
  13. 26598
  14. 4665
  15. 12958
  16. 12959
  17. 1194
  18. 21485
  19. 27734
  20. 15003
  21. 18904
  22. 16209
  23. 12296
  24. 13230
  25. 25407
  26. 20545
  27. 15605
In [19]:
## calculate gene list module score for each cell
Ethylene_module_score <- as.numeric(apply(apply(as.matrix(rc.integrated@assays$SCT@data)[known.good.markers$Gene.ID,], 1, zscore), 1, mean))

Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 19.4 GiB”
In [20]:
## Add to meta data
rc.integrated$Ethylene_module_score <- Ethylene_module_score
In [27]:
FeaturePlot(rc.integrated, reduction = "umap", feature="Ethylene_module_score", order=TRUE, pt.size=1,cols=c("lightgray","red"))
In [31]:
table(rc.integrated$celltype.fine.anno)

Stem Cell Niche    Sclerenchyma        Root Cap    Atrichoblast     Trichoblast 
           7746            8493            4724            4675            5401 
         Cortex      Endodermis          Phloem       Pericycle       Exodermis 
          17257            4362            2455            5088           10326 
     Protoxylem       Metaxylem 
           7928            1288
In [32]:
## Prepare dataframe for heatmap
df <- data.frame(celltype=rc.integrated$celltype.fine.anno, ethylene=rc.integrated$Ethylene_module_score)
In [39]:
df %>% group_by(celltype) %>%
    dplyr::summarize(Ethylene_Module_Score = mean(ethylene, na.rm=TRUE))
A tibble: 12 × 2
celltypeEthylene_Module_Score
<fct><dbl>
Stem Cell Niche-0.05798506
Sclerenchyma -0.02912666
Root Cap 0.18970615
Atrichoblast -0.16937811
Trichoblast -0.13884373
Cortex 0.07530903
Endodermis -0.11585721
Phloem 0.02853841
Pericycle 0.04115489
Exodermis -0.03385258
Protoxylem 0.07034536
Metaxylem 0.04678408
In [29]:
colnames(rc.integrated@meta.data)
  1. 'orig.ident'
  2. 'nCount_RNA'
  3. 'nFeature_RNA'
  4. 'nCount_spliced_RNA'
  5. 'nFeature_spliced_RNA'
  6. 'nCount_unspliced_RNA'
  7. 'nFeature_unspliced_RNA'
  8. 'percent.mt'
  9. 'percent.cp'
  10. 'nCount_SCT'
  11. 'nFeature_SCT'
  12. 'pANN_0.25_0.15_2850'
  13. 'DF.classifications_0.25_0.15_2850'
  14. 'nCount_spliced_SCT'
  15. 'nFeature_spliced_SCT'
  16. 'nCount_unspliced_SCT'
  17. 'nFeature_unspliced_SCT'
  18. 'SCT_snn_res.0.5'
  19. 'seurat_clusters'
  20. 'pANN_0.25_0.15_2547'
  21. 'DF.classifications_0.25_0.15_2547'
  22. 'pANN_0.25_0.15_205'
  23. 'DF.classifications_0.25_0.15_205'
  24. 'celltype.plate.ID'
  25. 'timezone.ID'
  26. 'celltype.plate.cor'
  27. 'timezone.cor'
  28. 'celltype.plate.pvalue'
  29. 'timezone.pvalue'
  30. 'celltype.field.ID'
  31. 'celltype.field.cor'
  32. 'celltype.field.pvalue'
  33. 'celltype.plate.ID.P'
  34. 'timezone.ID.P'
  35. 'celltype.plate.cor.P'
  36. 'timezone.cor.P'
  37. 'celltype.plate.pvalue.P'
  38. 'timezone.pvalue.P'
  39. 'celltype.field.ID.P'
  40. 'celltype.field.cor.P'
  41. 'celltype.field.pvalue.P'
  42. 'SCT_snn_res.50'
  43. 'seucluster.celltype.plate.cor.P'
  44. 'seucluster.celltype.plate.pvalue.P'
  45. 'seucluster.celltype.plate.ID.P'
  46. 'seucluster.celltype.field.cor.P'
  47. 'seucluster.celltype.field.pvalue.P'
  48. 'seucluster.celltype.field.ID.P'
  49. 'seucluster.timezone.cor.P'
  50. 'seucluster.timezone.pvalue.P'
  51. 'seucluster.timezone.ID.P'
  52. 'seucluster.celltype.plate.cor'
  53. 'seucluster.celltype.plate.pvalue'
  54. 'seucluster.celltype.plate.ID'
  55. 'seucluster.celltype.field.cor'
  56. 'seucluster.celltype.field.pvalue'
  57. 'seucluster.celltype.field.ID'
  58. 'seucluster.timezone.cor'
  59. 'seucluster.timezone.pvalue'
  60. 'seucluster.timezone.ID'
  61. 'pANN_0.25_0.15_166'
  62. 'DF.classifications_0.25_0.15_166'
  63. 'pANN_0.25_0.15_75'
  64. 'DF.classifications_0.25_0.15_75'
  65. 'SCT_snn_res.5'
  66. 'pANN_0.25_0.15_93'
  67. 'DF.classifications_0.25_0.15_93'
  68. 'pANN_0.25_0.15_222'
  69. 'DF.classifications_0.25_0.15_222'
  70. 'pANN_0.25_0.15_371'
  71. 'DF.classifications_0.25_0.15_371'
  72. 'pANN_0.25_0.15_232'
  73. 'DF.classifications_0.25_0.15_232'
  74. 'pANN_0.25_0.15_174'
  75. 'DF.classifications_0.25_0.15_174'
  76. 'pANN_0.25_0.15_190'
  77. 'DF.classifications_0.25_0.15_190'
  78. 'pANN_0.25_0.15_170'
  79. 'DF.classifications_0.25_0.15_170'
  80. 'integrated_snn_res.0.5'
  81. 'DF.classifications'
  82. 'log_nCount_RNA'
  83. 'geno'
  84. 'integrated_snn_res.2'
  85. 'score.crude.anno'
  86. 'integrated_snn_res.300'
  87. 'score.anno'
  88. 'consensus.anno'
  89. 'new_ref'
  90. 'final.anno.P'
  91. 'final.anno.cor.P'
  92. 'final.anno.pvalue.P'
  93. 'celltype.anno'
  94. 'time.anno'
  95. 'celltype.fine.anno'
  96. 'integrated_snn_res.500'
  97. 'Ethylene_module_score'
In [12]:
# Find included samples
table(rc.integrated$orig.ident)

sc_108 sc_109 sc_115 sc_116 sc_192 sc_193 sc_194 sc_195 sc_196   sc_7    tz1 
  4995   4499   3381   5188   6650   5296   4608   4802   4549   3042  15979 
   tz2 
 16754
In [5]:
# Check what variables are included in the metadata
head(rc.integrated@meta.data)
A data.frame: 6 × 96
orig.identnCount_RNAnFeature_RNAnCount_spliced_RNAnFeature_spliced_RNAnCount_unspliced_RNAnFeature_unspliced_RNApercent.mtpercent.cpnCount_SCT⋯score.annoconsensus.annonew_reffinal.anno.Pfinal.anno.cor.Pfinal.anno.pvalue.Pcelltype.annotime.annocelltype.fine.annointegrated_snn_res.500
<chr><dbl><int><dbl><int><dbl><int><dbl><dbl><dbl>⋯<fct><fct><chr><chr><dbl><dbl><fct><chr><fct><fct>
AAACCCAAGAATACAC_1tz2 75332277 72352203 2981410.464622330.013274924 9654⋯Endodermis Endodermis Maturation_Endodermis Maturation_Endodermis 0.93105451.228848e-88Endodermis MaturationEndodermis 341
AAACCCAAGAGTTCGG_1tz2 64361952 60691868 3671310.528278430.000000000 9882⋯SclerenchymaSclerenchymaMaturation_SclerenchymaMaturation_Sclerenchyma0.83734727.925984e-54SclerenchymaMaturationSclerenchyma6
AAACCCAAGGTAGCAC_1tz2 61782392 59022268 2761820.728391070.016186468 9468⋯Exodermis Cortex Maturation_Cortex Elongation_Exodermis 0.47285771.548340e-12Exodermis ElongationExodermis 2067
AAACCCAAGTGATGGC_1tz228987489527735468412524220.420878320.00000000010750⋯Exodermis Exodermis Elongation_Exodermis Elongation_Exodermis 0.86399426.486327e-61Exodermis ElongationExodermis 750
AAACCCAAGTGGAAGA_1tz241734647539910618118246120.064695450.00239612810108⋯Exodermis Exodermis Maturation_Exodermis Maturation_Exodermis 0.88925093.514207e-69Exodermis MaturationExodermis 2602
AAACCCAAGTTGCTCA_1tz2 95793099 90492918 5302900.480217140.010439503 9931⋯Trichoblast Trichoblast Maturation_Trichoblast Maturation_Trichoblast 0.92129413.704953e-83Trichoblast MaturationTrichoblast 350
In [13]:
table(rc.integrated$celltype.fine.anno)

Stem Cell Niche    Sclerenchyma        Root Cap    Atrichoblast     Trichoblast 
           7746            8493            4724            4675            5401 
         Cortex      Endodermis          Phloem       Pericycle       Exodermis 
          17257            4362            2455            5088           10326 
     Protoxylem       Metaxylem 
           7928            1288
In [7]:
library(ggplot2)
In [8]:
# Make cell type annotation figure
order <- c("Root Cap","Stem Cell Niche", "Atrichoblast", "Trichoblast", "Exodermis","Sclerenchyma","Cortex", "Endodermis","Pericycle", "Phloem", "Protoxylem","Metaxylem")
palette <- c("#6bff00", "#360167","#59adf6", "#Ff4061", "#F8f38d", "#32CD32", "#BF40BF", "#Ffb480","#08cad1", "#9d94ff", "#B5651D", "#834200")
rc.integrated$celltype.fine.anno <- factor(rc.integrated$celltype.fine.anno , levels = order[sort(match(unique(rc.integrated$celltype.fine.anno),order))]) 
color <- palette[sort(match(unique(rc.integrated$celltype.fine.anno),order))]
options(repr.plot.width=36, repr.plot.height=36)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.fine.anno", cols= palette, pt.size = 2.0)+ggtitle("Celltype Annotation UMAP")
In [ ]:
options(repr.plot.width=36, repr.plot.height=36)
DimPlot(rc.integrated, reduction = "umap", group.by = "time.anno", pt.size = 2.0)+ggtitle("Stage Annotation UMAP")

Correlation analysis¶

In [55]:
names(table(rc.integrated$orig.ident))
  1. 'sc_108'
  2. 'sc_109'
  3. 'sc_115'
  4. 'sc_116'
  5. 'sc_192'
  6. 'sc_193'
  7. 'sc_194'
  8. 'sc_195'
  9. 'sc_196'
  10. 'sc_7'
  11. 'tz1'
  12. 'tz2'
In [56]:
dat <- c()
for (i in names(table(rc.integrated$orig.ident))){
  dat <- cbind(dat, rowSums(as.matrix(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident==i)]))/ncol(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident==i)]))
}

Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.2 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.1 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.3 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.6 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.3 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.1 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.2 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.1 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 3.9 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 4.1 GiB”
In [58]:
colnames(dat) <- names(table(rc.integrated$orig.ident))
head(dat)
A matrix: 6 × 12 of type dbl
sc_108sc_109sc_115sc_116sc_192sc_193sc_194sc_195sc_196sc_7tz1tz2
LOC-Os01g010100.015581480.0149621010.0200255570.0174824940.0088795250.012545130.014003550.0097202060.0121268900.013349480.0022090490.02028568
LOC-Os01g010190.014884920.0135085430.0044600230.0127506200.0124041950.013150370.013226120.0138196720.0123825830.019600330.0035611580.01202689
LOC-Os01g010300.013070500.0072557420.0051949850.0098440390.0112414530.011147360.014679230.0102867350.0085446220.013766480.0024153090.02049023
LOC-Os01g010400.141101640.1081483920.1284486150.1751091700.1556812060.173011220.244428150.2102412370.2262946100.192420560.0445887510.23713867
LOC-Os01g010500.035993550.0422820510.0430415470.0543655460.0440278880.052318470.077061810.0705559270.0307684920.072938330.0047004830.01852074
LOC-Os01g010600.295548210.3826961410.1733443590.1592618560.3544814630.362472960.341335540.2833168630.2135188100.316248550.1128618430.45837298
In [59]:
options(repr.plot.width=12, repr.plot.height=12)
##The agglomeration method to be used. This should be (an unambiguous abbreviation of) one of "ward.D", "ward.D2", "single", "complete", "average" (= UPGMA), "mcquitty" (= WPGMA), "median" (= WPGMC) or "centroid" (= UPGMC).
##The distance measure to be used. This must be one of "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski". Any unambiguous substring can be given.
gplots::heatmap.2(cor(dat,method = "spearman"),col = brewer.pal(11,"YlOrRd"),  
          hclustfun = function(x) hclust(x, method = "ward.D2"), dendrogram="none",
          distfun = function(x) dist(x, method = "euclidean"), key.title = "Spearman correlation", trace = "none",cexCol = 1, cexRow = 1, margin=c(14,14),
                  Colv=FALSE, Rowv = FALSE)

Warning message in brewer.pal(11, "YlOrRd"):
“n too large, allowed maximum for palette YlOrRd is 9
Returning the palette you asked for with that many colors
”
In [15]:
# Find cell type markers to check the overlap the genes 
DefaultAssay(rc.integrated) <- "SCT"
Idents(rc.integrated) <- rc.integrated$celltype.fine.anno
rc.markers <- FindAllMarkers(rc.integrated, only.pos = TRUE, test.use = "wilcox",min.pct = 0.25, logfc.threshold = 0.25)
head(rc.markers, n =7)

Calculating cluster Root Cap

Calculating cluster Stem Cell Niche

Calculating cluster Atrichoblast

Calculating cluster Trichoblast

Calculating cluster Exodermis

Calculating cluster Sclerenchyma

Calculating cluster Cortex

Calculating cluster Endodermis

Calculating cluster Pericycle

Calculating cluster Phloem

Calculating cluster Protoxylem

Calculating cluster Metaxylem
A data.frame: 7 × 7
p_valavg_logFCpct.1pct.2p_val_adjclustergene
<dbl><dbl><dbl><dbl><dbl><fct><chr>
LOC-Os11g3797002.3573650.5540.1060Root CapLOC-Os11g37970
LOC-Os01g7320002.3300380.5780.1370Root CapLOC-Os01g73200
LOC-Os01g2235202.1921860.7150.1640Root CapLOC-Os01g22352
LOC-Os03g1885002.0317490.8320.3990Root CapLOC-Os03g18850
LOC-Os07g0346702.0302730.3180.0510Root CapLOC-Os07g03467
LOC-Os10g2179002.0171820.6810.2480Root CapLOC-Os10g21790
LOC-Os01g2237001.9875680.7100.2300Root CapLOC-Os01g22370
In [16]:
# Save all the filtered DE genes
write.csv(rc.markers, file = "./1-Root_Atlas/Cell_type_markers.csv", row.names = TRUE)
In [39]:
# Check marker gene expression-Root cap
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=18, repr.plot.height=6)
rootcap.marker <- c("LOC-Os03g18130","LOC-Os03g14300")
FeaturePlot(rc.integrated, features = rootcap.marker,
            cols = c("grey", "#6bff00"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order=T,ncol = 3)
In [5]:
# Check marker gene expression-Stem cell niche/meristem
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=18, repr.plot.height=12)
stemcellniche.marker <- c("LOC-Os01g16650","LOC-Os02g56130","LOC-Os03g51970","LOC-Os04g41900","LOC-Os05g44400")
FeaturePlot(rc.integrated, features = stemcellniche.marker,
            cols = c("grey", "#360167"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order=T, ncol = 3)
In [62]:
# Check marker gene expression-Epidermis
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=14, repr.plot.height=6)
epidermis.marker <- c("LOC-Os03g19990")
FeaturePlot(rc.integrated, features = epidermis.marker,
            cols = c("grey", "#c6b7fe"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order=T, ncol = 3)
In [26]:
# Check marker gene expression-Atrichoblast
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=18, repr.plot.height=6)
atrichoblast.marker <- c("LOC-Os01g14650","LOC-Os01g50820","LOC-Os01g64840")
FeaturePlot(rc.integrated, features = atrichoblast.marker,
            cols = c("grey", "#59adf6"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order=T, ncol = 3)
In [27]:
# Check marker gene expression-Trichoblast
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=18, repr.plot.height=24)
trichoblast.marker <- c("LOC-Os12g05380","LOC-Os05g38770","LOC-Os10g03400",
                        "LOC-Os04g45290","LOC-Os06g48050","LOC-Os01g03530",
                        "LOC-Os01g11750","LOC-Os03g04210","LOC-Os03g19330",
                        "LOC-Os10g42750","LOC-Os06g08500")
FeaturePlot(rc.integrated, features = trichoblast.marker,
            cols = c("grey", "#Ff4061"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order=T, ncol = 3)
In [71]:
# Check marker gene expression-Trichoblast
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, features = "LOC-Os05g38770",
            cols = c("grey", "#Ff4061"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order =T)
In [20]:
# Check marker gene expression-plot both trichoblast and trichoblast
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=32, repr.plot.height=8)
FeaturePlot(rc.integrated, features = c("LOC-Os05g38770", "LOC-Os01g64840"),
            cols = c("grey", "#Ff4061","#59adf6"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order =T, blend = TRUE)
In [29]:
# Check marker gene expression-Exodermis
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=18, repr.plot.height=12)
exodermis.marker <- c("LOC-Os03g37411","LOC-Os04g37980","LOC-Os03g02460",
                      "LOC-Os06g17260")
FeaturePlot(rc.integrated, features = exodermis.marker,
            cols = c("grey", "#F8f38d"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order=T, ncol = 3)
In [70]:
# Check marker gene expression-Trichoblast
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, features = "LOC-Os03g37411",
            cols = c("grey", "#F8f38d"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order =T)
In [31]:
# Check marker gene expression-Sclerenchyma
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=18, repr.plot.height=6)
sclerenchyma.marker <- c("LOC-Os08g02300","LOC-Os05g46610","LOC-Os08g05520")
FeaturePlot(rc.integrated, features = sclerenchyma.marker,
            cols = c("grey", "#32CD32"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order=T, ncol = 3)
In [69]:
# Check marker gene expression-Trichoblast
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, features = "LOC-Os08g02300",
            cols = c("grey", "#32CD32"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order = T)
In [32]:
# Check marker gene expression-Cortex
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=18, repr.plot.height=12)
cortex.marker <- c("LOC-Os03g04310","LOC-Os06g30730","LOC-Os05g33080",
                   "LOC-Os01g19220","LOC-Os04g46810")
FeaturePlot(rc.integrated, features = cortex.marker,
            cols = c("grey", "#BF40BF"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order=T, ncol = 3)
In [72]:
# Check marker gene expression-Trichoblast
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, features = "LOC-Os03g04310",
            cols = c("grey", "#BF40BF"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order = T)
In [34]:
# Check marker gene expression-Endodermis
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=18, repr.plot.height=12)
endodermis.marker <- c("LOC-Os01g15810","LOC-Os03g18640","LOC-Os10g06680",
                       "LOC-Os01g67390")
FeaturePlot(rc.integrated, features = endodermis.marker,
            cols = c("grey", "#ffb480"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order=T, ncol = 3)
In [73]:
# Check marker gene expression-Trichoblast
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, features = "LOC-Os01g15810",
            cols = c("grey", "#ffb480"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order = T)
In [35]:
# Check marker gene expression-Pericycle
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=18, repr.plot.height=12)
pericycle.marker <- c("LOC-Os02g56510","LOC-Os01g19170","LOC-Os01g58910",
                      "LOC-Os02g52930","LOC-Os07g44060")
FeaturePlot(rc.integrated, features = pericycle.marker,
            cols = c("grey", "#08cad1"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order=T, ncol = 3)
In [36]:
# Check marker gene expression-Phloem
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=18, repr.plot.height=6)
phloem.marker <- c("LOC-Os08g04400","LOC-Os06g45410","LOC-Os01g52480")
FeaturePlot(rc.integrated, features = phloem.marker,
            cols = c("grey", "#9d94ff"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order=T, ncol = 3)
In [74]:
# Check marker gene expression-Trichoblast
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, features = "LOC-Os06g45410",
            cols = c("grey", "#9d94ff"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order = T)
In [37]:
# Check marker gene expression-Xylem
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=18, repr.plot.height=12)
xylem.marker <- c("LOC-Os01g67090","LOC-Os01g48130","LOC-Os09g25490",
                  "LOC-Os10g32980","LOC-Os01g54620")
FeaturePlot(rc.integrated, features = xylem.marker,
            cols = c("grey", "#B7410E"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order=T, ncol = 3)
In [38]:
# Check marker gene expression-Metaxylem
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=18, repr.plot.height=6)
metaxylem.marker <- c("LOC-Os01g73980","LOC-Os07g44450","LOC-Os08g38170")
FeaturePlot(rc.integrated, features = metaxylem.marker,
            cols = c("grey", "#834200"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order=T, ncol = 3)
In [48]:
# Try the dotplot for the marker genes
all.marker <- read.csv("./1-Root_Atlas//Marker_genes_list_for-dotplot.csv")
all.marker$Gene.ID
  1. 'LOC-Os03g18130'
  2. 'LOC-Os03g14300'
  3. 'LOC-Os01g16650'
  4. 'LOC-Os02g56130'
  5. 'LOC-Os03g51970'
  6. 'LOC-Os04g41900'
  7. 'LOC-Os05g44400'
  8. 'LOC-Os03g19990'
  9. 'LOC-Os01g50820'
  10. 'LOC-Os01g64840'
  11. 'LOC-Os01g14650'
  12. 'LOC-Os12g05380'
  13. 'LOC-Os05g38770'
  14. 'LOC-Os10g03400'
  15. 'LOC-Os04g45290'
  16. 'LOC-Os06g48050'
  17. 'LOC-Os01g03530'
  18. 'LOC-Os01g11750'
  19. 'LOC-Os03g04210'
  20. 'LOC-Os03g19330'
  21. 'LOC-Os10g42750'
  22. 'LOC-Os06g08500'
  23. 'LOC-Os03g37411'
  24. 'LOC-Os04g37980'
  25. 'LOC-Os03g02460'
  26. 'LOC-Os06g17260'
  27. 'LOC-Os08g02300'
  28. 'LOC-Os05g46610'
  29. 'LOC-Os08g05520'
  30. 'LOC-Os03g04310'
  31. 'LOC-Os06g30730'
  32. 'LOC-Os05g33080'
  33. 'LOC-Os01g19220'
  34. 'LOC-Os04g46810'
  35. 'LOC-Os01g15810'
  36. 'LOC-Os03g18640'
  37. 'LOC-Os10g06680'
  38. 'LOC-Os01g67390'
  39. 'LOC-Os02g56510'
  40. 'LOC-Os01g19170'
  41. 'LOC-Os01g58910'
  42. 'LOC-Os02g52930'
  43. 'LOC-Os07g44060'
  44. 'LOC-Os08g04400'
  45. 'LOC-Os06g45410'
  46. 'LOC-Os01g52480'
  47. 'LOC-Os01g67090'
  48. 'LOC-Os01g48130'
  49. 'LOC-Os09g25490'
  50. 'LOC-Os10g32980'
  51. 'LOC-Os01g54620'
  52. 'LOC-Os01g73980'
  53. 'LOC-Os07g44450'
  54. 'LOC-Os08g38170'
In [51]:
options(repr.plot.width=18, repr.plot.height=6)
Idents(rc.integrated) <- rc.integrated$celltype.fine.anno
DotPlot(rc.integrated,assay = "SCT",
  features = all.marker$Gene.ID,
  cols = c("lightgrey", "#360167"),
  col.min = -2.5,
  col.max = 2.5,
  dot.min = 0,
  dot.scale = 8,
  group.by = NULL,
  split.by = NULL,
  scale = TRUE,
  scale.by = "radius",
  scale.min = NA,
  scale.max = NA
)
In [ ]:

Check the expression of specific genes¶

In [4]:
# Check marker gene expression-Stem cell niche/meristem
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=18, repr.plot.height=12)
marker <- c("LOC-Os02g47640","LOC-Os02g02690", "LOC-Os11g26910",
            "LOC-Os09g10260")
FeaturePlot(rc.integrated, features = marker,
            cols = c("grey", "#360167"), label=F, repel=F, pt.size = 0.4, 
            label.size = 1.5, order=T, ncol = 3)
In [ ]: